use "HSSC replication data\Generated data\model_results_saved.dta"model_results_saved.dta", clear

xtile lock10=RTstringencyindex, nq(100)
collapse (semean) sePharmincQ1=PharmincQ1 sePharmincQ2=PharmincQ2 sePharmincQ3=PharmincQ3 sePharmincQ4=PharmincQ4 sePharmincQ5=PharmincQ5 ///
sePjobincQ1=PjobincQ1 sePjobincQ2=PjobincQ2 sePjobincQ3=PjobincQ3 sePjobincQ4=PjobincQ4 sePjobincQ5=PjobincQ5 ///
sePincincQ1=PincincQ1 sePincincQ2=PincincQ2 sePincincQ3=PincincQ3 sePincincQ4=PincincQ4 sePincincQ5=PincincQ5 ///
(mean)  Pharm* Pjob* Paffect* Ptemp* Pinc* Phours* [aw=WGT], by(lock10)

reshape long sePjobinc sePharminc Pharminc Pjobinc Paffectinc Ptempinc Pincinc Phoursinc, i(lock10 ) j(inc) string

cd "HSSC replication data\Findings"

gen upHarm=Pharminc+(sePharminc*1.96)
gen lowHarm=Pharminc-(sePharminc*1.96)
gen upJob=Pjobinc+(sePjobinc*1.96)
gen lowJob=Pjobinc-(sePjobinc*1.96)

foreach x in Pjobinc Paffectinc Ptempinc Pincinc Phoursinc upHarm lowHarm upJob lowJob {
replace `x'=`x'*100
}

twoway lfitci Pharminc lock10 if inc=="Q1", clpattern(solid) clcolor(blue*.7) acolor(blue*.7) fintensity(inten10)  alpattern(dot) clwidth(vthick) ///
|| lfitci Pharminc lock10 if inc=="Q2", clpattern(dash) clcolor(blue*.5)  fintensity(inten10) alpattern(dot) clwidth(medthick) ///
|| lfitci Pharminc lock10 if inc=="Q3", clpattern(shortdash) clcolor(green*.7) fintensity(inten10) alpattern(dot) clwidth(medium) ///
|| lfitci Pharminc lock10 if inc=="Q4", clpattern(dash_dot) clcolor(red*.5) fintensity(inten10) alpattern(dot) clwidth(medthick) ///
|| lfitci Pharminc lock10 if inc=="Q5", clpattern(longdash) clcolor(red*.7) acolor(red*.7) fintensity(inten10)  alpattern(dot) clwidth(vthick) ///
|| scatter Pharminc lock10 if inc=="Q1", msymbol(dh) mcolor(blue*.7) ///
|| scatter Pharminc lock10 if inc=="Q5", msymbol(oh) mcolor(red*.7)  ///
plotregion(color(white)) graphregion(color(white)) ///
saving(Int1) ///
title("A. Harm index") ///
legend(order(2 "Bottom quintile" 4 "Quintile 2" 6 "Quintile 3" 8 "Quintile 4" 10 "Top quintile")    region(style(none))  size(small) rows(2) span)  ///
ytitle("Harm index") xtitle("Stringency in centiles") xlabel(0(20)100,  grid gmax labgap(tiny))

twoway lfitci Pjobinc lock10 if inc=="Q1", clpattern(solid) clcolor(blue*.7) acolor(blue*.7) fintensity(inten10)  alpattern(dot) clwidth(vthick) ///
|| lfitci Pjobinc lock10 if inc=="Q2", clpattern(dash) clcolor(blue*.5)  fintensity(inten10) alpattern(dot) clwidth(medthick) ///
|| lfitci Pjobinc lock10 if inc=="Q3", clpattern(shortdash) clcolor(green*.7) fintensity(inten10) alpattern(dot) clwidth(medium) ///
|| lfitci Pjobinc lock10 if inc=="Q4", clpattern(dash_dot) clcolor(red*.5) fintensity(inten10) alpattern(dot) clwidth(medthick) ///
|| lfitci Pjobinc lock10 if inc=="Q5", clpattern(longdash) clcolor(red*.7) acolor(red*.7) fintensity(inten10)  alpattern(dot) clwidth(vthick) ///
|| scatter Pjobinc lock10 if inc=="Q1", msymbol(dh) mcolor(blue*.7) ///
|| scatter Pjobinc lock10 if inc=="Q5", msymbol(oh) mcolor(red*.7)  ///
plotregion(color(white)) graphregion(color(white)) ///
saving(Int2) ///
title("B. Job loss") ///
legend(order(2 "Bottom quintile" 4 "Quintile 2" 6 "Quintile 3" 8 "Quintile 4" 10 "Top quintile")    region(style(none))  size(small) rows(2) span)  ///
ytitle("% lost job") xtitle("Stringency in centiles") xlabel(0(20)100,  grid gmax labgap(tiny))


twoway lfitci Paffectinc lock10 if inc=="Q1", clpattern(solid) clcolor(blue*.7) acolor(blue*.7) fintensity(inten10)  alpattern(dot) clwidth(vthick) ///
|| lfitci Paffectinc lock10 if inc=="Q2", clpattern(dash) clcolor(blue*.5)  fintensity(inten10) alpattern(dot) clwidth(medthick) ///
|| lfitci Paffectinc lock10 if inc=="Q3", clpattern(shortdash) clcolor(green*.7) fintensity(inten10) alpattern(dot) clwidth(medium) ///
|| lfitci Paffectinc lock10 if inc=="Q4", clpattern(dash_dot) clcolor(red*.5) fintensity(inten10) alpattern(dot) clwidth(medthick) ///
|| lfitci Paffectinc lock10 if inc=="Q5", clpattern(longdash) clcolor(red*.7) acolor(red*.7) fintensity(inten10)  alpattern(dot) clwidth(vthick) ///
|| scatter Paffectinc lock10 if inc=="Q1", msymbol(dh) mcolor(blue*.7) ///
|| scatter Paffectinc lock10 if inc=="Q5", msymbol(oh) mcolor(red*.7)  ///
plotregion(color(white)) graphregion(color(white)) ///
saving(Int3) ///
title("C. Affected a lot") ///
legend(order(2 "Bottom quintile" 4 "Quintile 2" 6 "Quintile 3" 8 "Quintile 4" 10 "Top quintile")    region(style(none))  size(small) rows(2) span)  ///
ytitle("% affected a lot") xtitle("Stringency in centiles") xlabel(0(20)100,  grid gmax labgap(tiny))

twoway lfitci Ptempinc lock10 if inc=="Q1", clpattern(solid) clcolor(blue*.7) acolor(blue*.7) fintensity(inten10)  alpattern(dot) clwidth(vthick) ///
|| lfitci Ptempinc lock10 if inc=="Q2", clpattern(dash) clcolor(blue*.5)  fintensity(inten10) alpattern(dot) clwidth(medthick) ///
|| lfitci Ptempinc lock10 if inc=="Q3", clpattern(shortdash) clcolor(green*.7) fintensity(inten10) alpattern(dot) clwidth(medium) ///
|| lfitci Ptempinc lock10 if inc=="Q4", clpattern(dash_dot) clcolor(red*.5) fintensity(inten10) alpattern(dot) clwidth(medthick) ///
|| lfitci Ptempinc lock10 if inc=="Q5", clpattern(longdash) clcolor(red*.7) acolor(red*.7) fintensity(inten10)  alpattern(dot) clwidth(vthick) ///
|| scatter Ptempinc lock10 if inc=="Q1", msymbol(dh) mcolor(blue*.7) ///
|| scatter Ptempinc lock10 if inc=="Q5", msymbol(oh) mcolor(red*.7)  ///
plotregion(color(white)) graphregion(color(white)) ///
saving(Int4) ///
title("D. Temporarily laid off") ///
legend(order(2 "Bottom quintile" 4 "Quintile 2" 6 "Quintile 3" 8 "Quintile 4" 10 "Top quintile")    region(style(none))  size(small) rows(2) span)  ///
ytitle("% temp laid off") xtitle("Stringency in centiles") xlabel(0(20)100,  grid gmax labgap(tiny))


twoway lfitci Pincinc lock10 if inc=="Q1", clpattern(solid) clcolor(blue*.7) acolor(blue*.7) fintensity(inten10)  alpattern(dot) clwidth(vthick) ///
|| lfitci Pincinc lock10 if inc=="Q2", clpattern(dash) clcolor(blue*.5)  fintensity(inten10) alpattern(dot) clwidth(medthick) ///
|| lfitci Pincinc lock10 if inc=="Q3", clpattern(shortdash) clcolor(green*.7) fintensity(inten10) alpattern(dot) clwidth(medium) ///
|| lfitci Pincinc lock10 if inc=="Q4", clpattern(dash_dot) clcolor(red*.5) fintensity(inten10) alpattern(dot) clwidth(medthick) ///
|| lfitci Pincinc lock10 if inc=="Q5", clpattern(longdash) clcolor(red*.7) acolor(red*.7) fintensity(inten10)  alpattern(dot) clwidth(vthick) ///
|| scatter Pincinc lock10 if inc=="Q1", msymbol(dh) mcolor(blue*.7) ///
|| scatter Pincinc lock10 if inc=="Q5", msymbol(oh) mcolor(red*.7)  ///
plotregion(color(white)) graphregion(color(white)) ///
saving(Int5) ///
title("D. Income loss") ///
legend(order(2 "Bottom quintile" 4 "Quintile 2" 6 "Quintile 3" 8 "Quintile 4" 10 "Top quintile")    region(style(none))  size(small) rows(2) span)  ///
ytitle("% lost income") xtitle("Stringency in centiles") xlabel(0(20)100,  grid gmax labgap(tiny))

grc1leg Int1.gph Int2.gph Int3.gph Int4.gph     ///
, imargin(.5 .5 .5 .5) row(2) col(2) legendfrom(Int1.gph) iscale(*.65) plotregion(color(white)) graphregion(color(white))  
graph export "HSSC replication data\Findings\FIG_Interaction_Models.png", replace

